#include <assert.h>
#include <stdio.h>

#include "../src/mstep.c"

static void
fill (double*           x,
      double          (*func) (void),
      unsigned int      n)
{
  for (unsigned int i = 0; i < n; ++i)
    {
      x[i] = func ();
    }
}

static void
fill_r (Rating*         x,
        unsigned int    n_workers,
        unsigned int    n_labels,
        unsigned int    n_images,
        unsigned int    n)
{
  for (unsigned int i = 0; i < n; ++i)
    {
      x[i].worker = (lrand48 () >> 16) % n_workers;
      x[i].image = (lrand48 () >> 16) % n_images;
      x[i].rating = (lrand48 () >> 16) % n_labels;
    }
}

static void
test_to_from_identity_once (void)
{
  /* (alpha, logbeta, gamma) = 
     from_x (alpha, logbeta, gamma, to_x (alpha, logbeta, gamma)) */

  unsigned int n_workers = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_labels = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_images = ((lrand48 () >> 16) % 5) + 1;

  double alpha[n_workers * n_labels * n_labels];
  double logbeta[n_images];
  double gamma[n_labels * n_labels];
  double newalpha[n_workers * n_labels * n_labels];
  double newlogbeta[n_images];
  double newgamma[n_labels * n_labels];

  fill (alpha, drand48, n_workers * n_labels * n_labels);
  fill (newalpha, drand48, n_workers * n_labels * n_labels);
  fill (logbeta, drand48, n_images);
  fill (newlogbeta, drand48, n_images);
  fill (gamma, drand48, n_labels * n_labels);
  fill (newgamma, drand48, n_labels * n_labels);

  size_t n_dims = 
    n_workers * n_labels * n_labels + n_images + n_labels * n_labels;
  gsl_vector* x = gsl_vector_alloc (n_dims);

  to_x (alpha, n_workers, n_labels, logbeta, n_images, gamma, x);
  from_x (newalpha, n_workers, n_labels, newlogbeta, n_images, newgamma, x);

  gsl_vector_free (x);

  assert (memcmp (alpha,
                  newalpha,
                  n_workers * n_labels * n_labels * sizeof (alpha[0])) == 0);

  assert (memcmp (logbeta, newlogbeta, n_images * sizeof (logbeta[0])) == 0);

  assert (memcmp (gamma, newgamma, n_labels * n_labels * sizeof (gamma[0])) == 0);
}

static void
test_to_from_inverse (void)
{
  for (unsigned int i = 0; i < 100; ++i)
    {
      test_to_from_identity_once ();
    }
}

static void
test_qfunc_f_once (void)
{
  /* qfunc (alpha, logbeta) = -qfunc_f (to_x (alpha, logbeta)) */

  unsigned int n_workers = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_labels = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_images = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_ratings = ((lrand48 () >> 16) % 5) + 1;

  double alpha[n_workers * n_labels * n_labels];
  double logbeta[n_images];
  double gamma[n_labels * n_labels];
  double pz[n_images * n_labels];
  double priorz[n_labels];
  Rating ratings[n_ratings];
  ParameterizedNormalDistribution prioralpha = parameterized_normal_distribution (1.0 + drand48 ());
  NormalDistribution priorlogbeta = normal_distribution (drand48 (), 1.0 + drand48 ());
  NormalDistribution priorgamma = normal_distribution (drand48 (), 1.0 + drand48 ());

  fill (alpha, drand48, n_workers * n_labels * n_labels);
  fill (logbeta, drand48, n_images);
  fill (gamma, drand48, n_labels * n_labels);
  fill (pz, drand48, n_images * n_labels);
  fill (priorz, drand48, n_labels);
  fill_r (ratings, n_workers, n_labels, n_images, n_ratings);

  double hq = 0;
  for (unsigned int k = 0; k < n_images * n_labels; ++k)
    {
      hq -= pz[k] * log (pz[k]);
    }

  QfuncParams qfunc_params = { .n_workers = n_workers,
                               .n_labels = n_labels,
                               .n_images = n_images,
                               .pz = pz,
                               .priorz = priorz,
                               .ratings = ratings,
                               .n_ratings = n_ratings,
                               .prioralpha = &prioralpha.base,
                               .priorlogbeta = &priorlogbeta.base,
                               .priorgamma = &priorgamma.base,
                               .symmetric = 0,
                               .hyperprior = 1,
                               .hq = hq };

  size_t n_dims = 
    n_workers * n_labels * n_labels + n_images + n_labels * n_labels;
  gsl_vector* x = gsl_vector_alloc (n_dims);

  to_x (alpha, n_workers, n_labels, logbeta, n_images, gamma, x);

  double q = qfunc (alpha, n_workers, n_labels, logbeta, n_images,
                    gamma, pz, ratings, n_ratings, 
                    &prioralpha.base, &priorlogbeta.base, &priorgamma.base,
                    0, hq);

  double qf = qfunc_f (x, &qfunc_params);

  assert (fabs (q + qf) <= 1e-8 * (1 + fabs (q) + fabs (qf)) ||
          (fprintf (stderr, "%g ?= %g\n", q, qf), 0));

  gsl_vector_free (x);
}

static void
test_qfunc_f (void)
{
  for (unsigned int i = 0; i < 100; ++i)
    {
      test_qfunc_f_once ();
    }
}

static void
test_qfunc_df_once (void)
{
  unsigned int n_workers = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_labels = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_images = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_ratings = ((lrand48 () >> 16) % 5) + 1;

  /* dqfunc (alpha, logbeta) = -qfunc_df (to_x (alpha, logbeta)) */

  double alpha[n_workers * n_labels * n_labels];
  double logbeta[n_images];
  double gamma[n_labels * n_labels];
  double dalpha[n_workers * n_labels * n_labels];
  double dlogbeta[n_images];
  double dgamma[n_labels * n_labels];
  double newdalpha[n_workers * n_labels * n_labels];
  double newdlogbeta[n_images];
  double newdgamma[n_labels * n_labels];
  double pz[n_images * n_labels];
  double priorz[n_labels];
  Rating ratings[n_ratings];
  ParameterizedNormalDistribution prioralpha = parameterized_normal_distribution (1.0 + drand48 ());
  NormalDistribution priorlogbeta = normal_distribution (drand48 (), 1.0 + drand48 ());
  NormalDistribution priorgamma = normal_distribution (drand48 (), 1.0 + drand48 ());

  fill (alpha, drand48, n_workers * n_labels * n_labels);
  fill (logbeta, drand48, n_images);
  fill (gamma, drand48, n_labels * n_labels);
  fill (dalpha, drand48, n_workers * n_labels * n_labels);
  fill (dlogbeta, drand48, n_images);
  fill (dgamma, drand48, n_labels * n_labels);
  fill (newdalpha, drand48, n_workers * n_labels * n_labels);
  fill (newdlogbeta, drand48, n_images);
  fill (newdgamma, drand48, n_labels * n_labels);
  fill (pz, drand48, n_images * n_labels);
  fill (priorz, drand48, n_labels);
  fill_r (ratings, n_workers, n_labels, n_images, n_ratings);

  double hq = 0;
  for (unsigned int k = 0; k < n_images * n_labels; ++k)
    {
      hq -= pz[k] * log (pz[k]);
    }

  QfuncParams qfunc_params = { .n_workers = n_workers,
                               .n_labels = n_labels,
                               .n_images = n_images,
                               .pz = pz,
                               .priorz = priorz,
                               .ratings = ratings,
                               .n_ratings = n_ratings,
                               .prioralpha = &prioralpha.base,
                               .priorlogbeta = &priorlogbeta.base,
                               .priorgamma = &priorgamma.base,
                               .symmetric = 0,
                               .hyperprior = 1,
                               .hq = hq };

  size_t n_dims = 
    n_workers * n_labels * n_labels + n_images + n_labels * n_labels;
  gsl_vector* x = gsl_vector_alloc (n_dims);
  gsl_vector* g = gsl_vector_alloc (n_dims);

  to_x (alpha, n_workers, n_labels, logbeta, n_images, gamma, x);

  qdqfunc (dalpha, dlogbeta, dgamma,
           alpha, n_workers, n_labels, logbeta, n_images, gamma,
           pz, ratings, n_ratings, 
           &prioralpha.base, &priorlogbeta.base, &priorgamma.base, 0, hq);

  qfunc_df (x, &qfunc_params, g);

  from_x (newdalpha, n_workers, n_labels, newdlogbeta, n_images, newdgamma, g);

  for (unsigned int i = 0; i < n_workers * n_labels * n_labels; ++i)
    {
      assert (   fabs (newdalpha[i] + dalpha[i]) 
              <= 1e-8 * (1 + fabs (newdalpha[i]) + fabs (dalpha[i])) ||
              (fprintf (stderr, "%u: %g ?= %g\n", i, newdalpha[i], dalpha[i]),
               0));
    }

  for (unsigned int i = 0; i < n_images; ++i)
    {
      assert (   fabs (newdlogbeta[i] + dlogbeta[i]) 
              <= 1e-8 * (1 + fabs (newdlogbeta[i]) + fabs (dlogbeta[i])) ||
              (fprintf (stderr, "%u: %g ?= %g\n", i, newdlogbeta[i], dlogbeta[i]),
               0));
    }

  for (unsigned int i = 0; i < n_labels * n_labels; ++i)
    {
      assert (   fabs (newdgamma[i] + dgamma[i]) 
              <= 1e-8 * (1 + fabs (newdgamma[i]) + fabs (dgamma[i])) ||
              (fprintf (stderr, "%u: %g ?= %g\n", i, newdgamma[i], dgamma[i]),
               0));
    }

  gsl_vector_free (g);
  gsl_vector_free (x);
}

static void
test_qfunc_df (void)
{
  for (unsigned int i = 0; i < 100; ++i)
    {
      test_qfunc_df_once ();
    }
}

static void
test_mstep_no_data_once (void)
{
  unsigned int n_workers = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_labels = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_images = ((lrand48 () >> 16) % 5) + 1;

  /* without data mstep should converge on prior maximum */

  double alpha[n_workers * n_labels * n_labels];
  double diagonal[n_workers * n_labels];
  double logbeta[n_images];
  double gamma[n_labels * n_labels];
  double gdiagonal[n_labels];
  double pz[n_images * n_labels];
  double priorz[n_labels];
  ParameterizedNormalDistribution prioralpha = parameterized_normal_distribution (1.0 + drand48 ());
  NormalDistribution priorlogbeta = normal_distribution (drand48 (), 1.0 + drand48 ());
  NormalDistribution priorgamma = normal_distribution (drand48 (), 1.0 + drand48 ());

  fill (alpha, drand48, n_workers * n_labels * n_labels);
  fill (logbeta, drand48, n_images);
  fill (gamma, drand48, n_labels * n_labels);
  fill (pz, drand48, n_images * n_labels);
  fill (priorz, drand48, n_labels);

  for (unsigned int j = 0; j < n_workers; ++j)
    {
      for (unsigned int k = 0; k < n_labels; ++k)
        {
          diagonal[j * n_labels + k] = 
            alpha[j * n_labels * n_labels + k * n_labels + k];
          gdiagonal[k] = gamma[k * n_labels + k];
        }
    }

  double hq = 0;
  for (unsigned int k = 0; k < n_images * n_labels; ++k)
    {
      hq -= pz[k] * log (pz[k]);
    }

  mstep (alpha, n_workers, n_labels, 
         logbeta, n_images, gamma,
         pz, priorz, 
         NULL, 0, 
         &prioralpha.base, &priorlogbeta.base, &priorgamma.base,
         0, 1,
         hq, 50, 1e-8);

  for (unsigned int i = 0, x = 0; i < n_workers * n_labels * n_labels; ++i)
    {
      if (i % n_labels == (i / n_labels) % n_labels)
        {
          assert (   fabs (alpha[i] - diagonal[x]) <= 
                     1e-8 * (1 + fabs (alpha[i]) + fabs (diagonal[x])) ||
                  (fprintf (stderr, "%u: %g ?= %g\n", i, alpha[i], diagonal[x]),
                   0));

          ++x;
        }
      else
        {
          assert (   fabs (alpha[i] - priorgamma.mean) <= 
                     1e-8 * (1 + fabs (alpha[i]) + fabs (priorgamma.mean)) ||
                  (fprintf (stderr, "%u: %g ?= %g\n", i, alpha[i], priorgamma.mean),
                   0));
        }
    }

  for (unsigned int i = 0; i < n_images; ++i)
    {
      assert (   fabs (logbeta[i] - priorlogbeta.mean) <= 
                 1e-8 * (1 + fabs (logbeta[i]) + fabs (priorlogbeta.mean)) ||
              (fprintf (stderr, "%u: %g ?= %g\n", i, logbeta[i], priorlogbeta.mean),
               0));
    }

  for (unsigned int i = 0, x = 0; i < n_labels * n_labels; ++i)
    {
      if (i % n_labels == (i / n_labels) % n_labels)
        {
          assert (   fabs (gamma[i] - gdiagonal[x]) <= 
                     1e-8 * (1 + fabs (gamma[i]) + fabs (gdiagonal[x])) ||
                  (fprintf (stderr, "%u: %g ?= %g\n", i, gamma[i], gdiagonal[x]),
                   0));

          ++x;
        }
      else
        {
          assert (   fabs (gamma[i] - priorgamma.mean) <= 
                     1e-8 * (1 + fabs (gamma[i]) + fabs (priorgamma.mean)) ||
                  (fprintf (stderr, "%u: %g ?= %g\n", i, gamma[i], priorgamma.mean),
                   0));
        }
    }
}

static void
test_mstep_no_data (void)
{
  for (unsigned int i = 0; i < 20; ++i)
    {
      test_mstep_no_data_once ();
    }
}

static void
test_mstep_symmetric_no_data_once (void)
{
  unsigned int n_workers = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_labels = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_images = ((lrand48 () >> 16) % 5) + 1;

  /* without data mstep should converge on prior maximum */

  double alpha[n_workers * n_labels * n_labels];
  double alphacopy[n_workers * n_labels * n_labels];
  double logbeta[n_images];
  double gamma[n_labels * n_labels];
  double gammacopy[n_labels * n_labels];
  double pz[n_images * n_labels];
  double priorz[n_labels];
  ParameterizedNormalDistribution prioralpha = parameterized_normal_distribution (1.0 + drand48 ());
  NormalDistribution priorlogbeta = normal_distribution (drand48 (), 1.0 + drand48 ());
  NormalDistribution priorgamma = normal_distribution (drand48 (), 1.0 + drand48 ());

  fill (alpha, drand48, n_workers * n_labels * n_labels);
  memcpy (alphacopy, alpha, n_workers * n_labels * n_labels * sizeof (alpha[0]));
  fill (logbeta, drand48, n_images);
  fill (gamma, drand48, n_labels * n_labels);
  memcpy (gammacopy, gamma, n_labels * n_labels * sizeof (gamma[0]));
  fill (pz, drand48, n_images * n_labels);
  fill (priorz, drand48, n_labels);

  double hq = 0;
  for (unsigned int k = 0; k < n_images * n_labels; ++k)
    {
      hq -= pz[k] * log (pz[k]);
    }

  mstep (alpha, n_workers, n_labels, 
         logbeta, n_images, gamma,
         pz, priorz, 
         NULL, 0, 
         &prioralpha.base, &priorlogbeta.base, &priorgamma.base, 
         1, 1,
         hq, 50, 1e-8);

  for (unsigned int i = 0; i < n_workers * n_labels * n_labels; ++i)
    {
      if (i % n_labels <= (i / n_labels) % n_labels)
        {
          assert (   fabs (alpha[i] - alphacopy[i]) <= 
                     1e-8 * (1 + fabs (alpha[i]) + fabs (alphacopy[i])) ||
                  (fprintf (stderr, "%u: %g ?= %g\n", i, alpha[i], alphacopy[i]),
                   0));
        }
      else
        {
          assert (   fabs (alpha[i] - priorgamma.mean) <= 
                     1e-8 * (1 + fabs (alpha[i]) + fabs (priorgamma.mean)) ||
                  (fprintf (stderr, "%u: %g ?= %g\n", i, alpha[i], priorgamma.mean),
                   0));
        }
    }

  for (unsigned int i = 0; i < n_images; ++i)
    {
      assert (   fabs (logbeta[i] - priorlogbeta.mean) <= 
                 1e-8 * (1 + fabs (logbeta[i]) + fabs (priorlogbeta.mean)) ||
              (fprintf (stderr, "%u: %g ?= %g\n", i, logbeta[i], priorlogbeta.mean),
               0));
    }

  for (unsigned int i = 0; i < n_labels * n_labels; ++i)
    {
      if (i % n_labels <= (i / n_labels) % n_labels)
        {
          assert (   fabs (gamma[i] - gammacopy[i]) <= 
                     1e-8 * (1 + fabs (gamma[i]) + fabs (gammacopy[i])) ||
                  (fprintf (stderr, "%u: %g ?= %g\n", i, gamma[i], gammacopy[i]),
                   0));
        }
      else
        {
          assert (   fabs (gamma[i] - priorgamma.mean) <= 
                     1e-8 * (1 + fabs (gamma[i]) + fabs (priorgamma.mean)) ||
                  (fprintf (stderr, "%u: %g ?= %g\n", i, gamma[i], priorgamma.mean),
                   0));
        }
    }
}

static void
test_mstep_symmetric_no_data (void)
{
  for (unsigned int i = 0; i < 20; ++i)
    {
      test_mstep_symmetric_no_data_once ();
    }
}

typedef struct _MstepTestDatum MstepTestDatum ;
struct _MstepTestDatum 
{ 
  struct
    {
      double    alpha[4 * 5 * 5];
      double    logbeta[3];
      double    gamma[5 * 5];
      double    pz[3 * 5];
      double    priorz[5];
      Rating    ratings[8];
      double    alphastddev;
      double    logbetamean;
      double    logbetastddev;
      double    gammamean;
      double    gammastddev;
    }                   inputs;

  struct
    {
      double    newq;
      struct
        {
          double    alpha[4 * 5 * 5];
          double    logbeta[3];
          double    gamma[5 * 5];
        }       new;
    }                   desired_output;
};

static void
test_mstep (void)
{
  MstepTestDatum data[] = 
    {{{{0, 0.951853, 0.99105, 0.737832, 0.58324, 0.258483, 0, 0.527755,0.595682, 0.023121, 0.365911, 0.424452, 0, 0.423352, 0.639832,0.93982, 0.2432, 0.172489, 0, 0.308446, 0.157189, 0.920438,0.841625, 0.0212579, 0, 0, 0.484874, 0.404654, 0.511396, 0.49086,0.53302, 0, 0.413604, 0.773564, 0.90762, 0.274537, 0.885849, 0,0.177882, 0.884498, 0.908626, 0.461397, 0.75453, 0, 0.244667,0.968806, 0.218196, 0.582041, 0.936221, 0, 0, 0.811617, 0.297758,0.740416, 0.914963, 0.326743, 0, 0.893104, 0.22902, 0.424103,0.793722, 0.4795, 0, 0.455456, 0.516484, 0.519185, 0.593651,0.277574, 0, 0.631985, 0.610559, 0.132255, 0.523044, 0.387318, 0,0, 0.641753, 0.914058, 0.941002, 0.451097, 0.830136, 0, 0.6163,0.200586, 0.536135, 0.503393, 0.723196, 0, 0.971566, 0.112031,0.709671, 0.243696, 0.516111, 0, 0.595548, 0.190485, 0.650045,0.238537, 0.963562, 0}, {0.579927, 0.51779, 0.715493}, {1.,0.576244, 0.938174, 0.603732, 0.774491, 0.125147, 1., 0.108038,0.987432, 0.573905, 0.589012, 0.604645, 1., 0.264235, 0.602339,0.476981, 0.894975, 0.0205393, 1., 0.0862282, 0.881433, 0.704489,0.370494, 0.847691, 1.}, {0.917871, 0.124563, 0.852704, 0.132198,0.341627, 0.186389, 0.248972, 0.357707, 0.21648, 0.0783507,0.261541, 0.783802, 0.627467, 0.473706, 0.997306}, {0.181463,0.150487, 0.578731, 0.976766,0.0952347}, {{3, 1, 1}, {3, 0, 3}, {2, 1, 1}, {0, 0, 4}, {1, 2,4}, {1, 2, 3}, {3, 0, 1}, {3, 0, 4}}, 0.269053, 0.874242,0.606272, 0.247543,0.351182}, {-130.225, {{0., 0.122056, 0.349565,0.105546, -0.213227, 0.335121, 0., 0.335121, 0.101749, 0.0244971,0.358838, 0.108041, 0., 0.0749929, -0.217469, 0.333236,0.197543, 0.333236, 0., 0.0643507, 0.382405, 0.261149, 0.382405,0.146445, 0., 0., 0.11235, 0.347417, 0.0177137, -0.0898168,0.35938, 0., 0.35938, -0.0162114, -0.0356967, 0.365576, 0.123349,0., -0.0619161, -0.1304, 0.353777, 0.232765, 0.353777, 0.,0.0429066, 0.414038, 0.310265, 0.414038, 0.0263395, 0., 0.,0.0600008, 0.342569, 0.0746111, -0.0556333, 0.337855, 0.,0.337855, 0.102111, 0.0684046, 0.356061, 0.0227751, 0.,0.0487214, -0.0594148, 0.334827, 0.15823, 0.334827, 0., 0.109214,0.375463, 0.23264, 0.375463, 0.119489, 0., 0., -0.0380758,0.371003, -0.0316901, -0.0847326, 0.34626, 0., 0.34626,0.0882512, 0.059835, 0.379644, -0.0563836,0., -0.0446868, -0.0910859, 0.344839, 0.138373, 0.344839, 0.,0.0966945, 0.396688, 0.172297, 0.396688, 0.0944614, 0.}, {1.6641,0.939379, 1.54596}, {1., 0.087559, 0.33919,0.0679053, -0.0649909, 0.332227, 1., 0.332227, 0.0918252,0.0571922, 0.349996, 0.0747946, 1., 0.0354069, -0.0769727,0.329625, 0.19015, 0.329625, 1., 0.0999495, 0.373644, 0.24453,0.373644, 0.115988, 1.}}}}};

  for (unsigned int i = 0; i < sizeof (data) / sizeof (data[0]); ++i)
    {
      ParameterizedNormalDistribution prior_alpha = 
        parameterized_normal_distribution (data[i].inputs.alphastddev);
      NormalDistribution prior_logbeta = 
        normal_distribution (data[i].inputs.logbetamean, 
                             data[i].inputs.logbetastddev);
      NormalDistribution prior_gamma = 
        normal_distribution (data[i].inputs.gammamean, 
                             data[i].inputs.gammastddev);
      double diagonal[4 * 5];
      double gdiagonal[5];

      for (unsigned int j = 0; j < 4; ++j)
        {
          for (unsigned int k = 0; k < 5; ++k)
            {
              diagonal[j * 5 + k] = drand48 ();
              data[i].inputs.alpha[j * 5 * 5 + k * 5 + k] = diagonal[j * 5 + k];
            }
        }

      for (unsigned int k = 0; k < 5; ++k)
        {
          gdiagonal[k] = drand48 ();
          data[i].inputs.gamma[k * 5 + k] = gdiagonal[k];
        }

      double hq = 0;
      for (unsigned int k = 0; k < 3 * 5; ++k)
        {
          hq -= data[i].inputs.pz[k] * log (data[i].inputs.pz[k]);
        }

      double q = mstep (data[i].inputs.alpha,
                        4,
                        5,
                        data[i].inputs.logbeta,
                        3,
                        data[i].inputs.gamma,
                        data[i].inputs.pz,
                        data[i].inputs.priorz,
                        data[i].inputs.ratings,
                        8,
                        &prior_alpha.base,
                        &prior_logbeta.base,
                        &prior_gamma.base,
                        0,
                        1,
                        hq,
                        1000,
                        1e-16);

      assert (fabs (q - data[i].desired_output.newq) <=
              1e-4 * (1 + fabs (q) + fabs (data[i].desired_output.newq)) ||
              ( fprintf (stderr, "%g ?= %g\n", q, data[i].desired_output.newq),
                0 ));
      
      for (unsigned int j = 0, y = 0; j < 4*5*5; ++j)
        {
          if (j % 5 == (j / 5) % 5)
            {
              assert (fabs (data[i].inputs.alpha[j] - diagonal[y]) <=
                      1e-8 * (1 + fabs (data[i].inputs.alpha[j]) + fabs (diagonal[y])) ||
                      ( fprintf (stderr, "%u %u: %g ?= %g\n", i, j, data[i].inputs.alpha[j], diagonal[y]),
                        0 ));

              ++y;
            }
          else
            {
              assert (fabs (data[i].inputs.alpha[j] - data[i].desired_output.new.alpha[j]) <=
                      5e-3 * (1 + fabs (data[i].inputs.alpha[j]) + fabs (data[i].desired_output.new.alpha[j])) ||
                      ( fprintf (stderr, "%u %u: %g ?= %g\n", i, j, data[i].inputs.alpha[j], data[i].desired_output.new.alpha[j]),
                        0 ));
            }
        }

      for (unsigned int j = 0; j < 3; ++j)
        {
          assert (fabs (data[i].inputs.logbeta[j] - data[i].desired_output.new.logbeta[j]) <=
                  5e-3 * (1 + fabs (data[i].inputs.logbeta[j]) + fabs (data[i].desired_output.new.logbeta[j])) ||
                  ( fprintf (stderr, "%u %u: %g ?= %g\n", i, j, data[i].inputs.logbeta[j], data[i].desired_output.new.logbeta[j]),
                    0 ));
        }

      for (unsigned int j = 0, y = 0; j < 5*5; ++j)
        {
          if (j % 5 == (j / 5) % 5)
            {
              assert (fabs (data[i].inputs.gamma[j] - gdiagonal[y]) <=
                      1e-8 * (1 + fabs (data[i].inputs.gamma[j]) + fabs (diagonal[y])) ||
                      ( fprintf (stderr, "%u %u: %g ?= %g\n", i, j, data[i].inputs.gamma[j], diagonal[y]),
                        0 ));

              ++y;
            }
          else
            {
              assert (fabs (data[i].inputs.gamma[j] - data[i].desired_output.new.gamma[j]) <=
                      5e-3 * (1 + fabs (data[i].inputs.gamma[j]) + fabs (data[i].desired_output.new.gamma[j])) ||
                      ( fprintf (stderr, "%u %u: %g ?= %g\n", i, j, data[i].inputs.gamma[j], data[i].desired_output.new.gamma[j]),
                        0 ));
            }
        }
    }
}

typedef struct _MstepNoHyperTestDatum MstepNoHyperTestDatum ;
struct _MstepNoHyperTestDatum 
{ 
  struct
    {
      double    alpha[4 * 5 * 5];
      double    logbeta[3];
      double    gamma[5 * 5];
      double    pz[3 * 5];
      double    priorz[5];
      Rating    ratings[8];
      double    alphastddev;
      double    logbetamean;
      double    logbetastddev;
      double    gammamean;
      double    gammastddev;
    }                   inputs;

  struct
    {
      double    newq;
      struct
        {
          double    alpha[4 * 5 * 5];
          double    logbeta[3];
        }       new;
    }                   desired_output;
};

static void
test_mstep_no_hyperprior (void)
{
  MstepNoHyperTestDatum data[] = 
    {{{{0, 0.677494, 0.925534, 0.853555, 0.0701316, 0.602037, 0, 0.0923837, 0.318849, 0.3969, 0.327083, 0.944772, 0, 0.50679, 0.59384, 0.96362, 0.641827, 0.282659, 0, 0.941948, 0.790273, 0.22151, 0.741002, 0.473761, 0, 0, 0.919146, 0.476158, 0.967302, 0.212698, 0.241652, 0, 0.550625, 0.113747, 0.142566, 0.639615, 0.458241, 0, 0.794898, 0.745666, 0.312532, 0.51347, 0.288108, 0, 0.151827, 0.348912, 0.871642, 0.00544943, 0.209879, 0, 0, 0.558639, 0.650133, 0.264447, 0.736118, 0.639493, 0, 0.173974, 0.297145, 0.52342, 0.397841, 0.623349, 0, 0.183399, 0.380855, 0.758226, 0.165108, 0.388501, 0, 0.635188, 0.445694, 0.651638, 0.100393, 0.483361, 0, 0, 0.0967821, 0.779996, 0.0949432, 0.273482, 0.538143, 0, 0.129863, 0.830496, 0.537364, 0.89865, 0.955889, 0, 0.53335, 0.0139437, 0.50081, 0.33254, 0.349952, 0, 0.633089, 0.742584, 0.167432, 0.961451, 0.997901, 0}, {0.29689, 0.515794, 0.861058}, {1., 0.51454, 0.200108, 0.735798, 0.766115, 0.241057, 1., 0.661964, 0.605934, 0.935619, 0.703693, 0.763314, 1., 0.650045, 0.402269, 0.68975, 0.262504, 0.317505, 1., 0.0523173, 0.0566604, 0.519921, 0.150073, 0.0908664, 1.}, {0.0587595, 0.223031, 0.63428, 0.229808, 0.54422, 0.022923, 0.898482, 0.463693, 0.303162, 0.360959, 0.292547, 0.528074, 0.599469, 0.597645, 0.642502}, {0.125805, 0.909719, 0.33514, 0.324997, 0.0734882}, {{3, 2, 1}, {3, 0, 0}, {1, 1, 0}, {0, 1, 3}, {3, 1, 0}, {3, 2, 4}, {2, 0, 1}, {3, 0, 3}}, 0.853059, 0.81522, 0.174924, 0.982622, 0.7943}, {-130.98, {{0., 0.519695, 0.210969, 0.698885, 0.76894, 0.419575, 0., 0.743629, -0.191638, 0.981604, 0.758587, 0.811614, 0., 0.102867, 0.503795, 0.724908, 0.3495, 0.395382, 0., 0.182495, 0.154981, 0.557169, 0.231763, -0.269311, 0., 0., 0.519735, 0.211051, 0.738859, 0.768962, -0.358743, 0., 0.732257, 0.684764, 0.974769, 0.139165, 0.812506, 0., 0.712643, 0.505473, 0.252452, 0.339588, 0.386361, 0., 0.168738, -0.292283, 0.556652, 0.23073, 0.181663, 0., 0., 0.432386, 0.225268, 0.743487, 0.773297, 0.328112, 0., 0.699212, 0.647828, 0.956305, 0.775009, 0.0673542, 0., 0.729233, 0.527808, 0.711991, -0.00530335, 0.366527, 0., 0.135481, 0.199107, -0.102243, 0.270935, 0.225092, 0., 0., 0.291827, 0.39567, 0.73079, 0.464736, -0.0694959, 0., 0.856779, 0.592253, 0.51547, 0.0751095, 0.397887, 0., 0.314446, 0.258711, 0.254674, 0.0559215, 0.559014, 0., -0.0178896, -0.091263, 0.203411, 0.508032, 0.0957806, 0.}, {0.827524, 0.88102, 0.790025}}}}};

  for (unsigned int i = 0; i < sizeof (data) / sizeof (data[0]); ++i)
    {
      ParameterizedNormalDistribution prior_alpha = 
        parameterized_normal_distribution (data[i].inputs.alphastddev);
      NormalDistribution prior_logbeta = 
        normal_distribution (data[i].inputs.logbetamean, 
                             data[i].inputs.logbetastddev);
      NormalDistribution prior_gamma = 
        normal_distribution (data[i].inputs.gammamean, 
                             data[i].inputs.gammastddev);
      double diagonal[4 * 5];
      double gammacopy[5 * 5];

      for (unsigned int j = 0; j < 4; ++j)
        {
          for (unsigned int k = 0; k < 5; ++k)
            {
              diagonal[j * 5 + k] = drand48 ();
              data[i].inputs.alpha[j * 5 * 5 + k * 5 + k] = diagonal[j * 5 + k];
            }
        }

      memcpy (gammacopy, data[i].inputs.gamma, 5 * 5 * sizeof (gammacopy[0]));

      double hq = 0;
      for (unsigned int k = 0; k < 3 * 5; ++k)
        {
          hq -= data[i].inputs.pz[k] * log (data[i].inputs.pz[k]);
        }

      double q = mstep (data[i].inputs.alpha,
                        4,
                        5,
                        data[i].inputs.logbeta,
                        3,
                        data[i].inputs.gamma,
                        data[i].inputs.pz,
                        data[i].inputs.priorz,
                        data[i].inputs.ratings,
                        8,
                        &prior_alpha.base,
                        &prior_logbeta.base,
                        &prior_gamma.base,
                        0,
                        0,
                        hq,
                        100,
                        1e-8);

      assert (fabs (q - data[i].desired_output.newq) <=
              1e-4 * (1 + fabs (q) + fabs (data[i].desired_output.newq)) ||
              ( fprintf (stderr, "%g ?= %g\n", q, data[i].desired_output.newq),
                0 ));
      
      for (unsigned int j = 0, y = 0; j < 4*5*5; ++j)
        {
          if (j % 5 == (j / 5) % 5)
            {
              assert (fabs (data[i].inputs.alpha[j] - diagonal[y]) <=
                      1e-8 * (1 + fabs (data[i].inputs.alpha[j]) + fabs (diagonal[y])) ||
                      ( fprintf (stderr, "%u %u: %g ?= %g\n", i, j, data[i].inputs.alpha[j], diagonal[y]),
                        0 ));

              ++y;
            }
          else
            {
              assert (fabs (data[i].inputs.alpha[j] - data[i].desired_output.new.alpha[j]) <=
                      5e-3 * (1 + fabs (data[i].inputs.alpha[j]) + fabs (data[i].desired_output.new.alpha[j])) ||
                      ( fprintf (stderr, "%u %u: %g ?= %g\n", i, j, data[i].inputs.alpha[j], data[i].desired_output.new.alpha[j]),
                        0 ));
            }
        }

      for (unsigned int j = 0; j < 3; ++j)
        {
          assert (fabs (data[i].inputs.logbeta[j] - data[i].desired_output.new.logbeta[j]) <=
                  5e-3 * (1 + fabs (data[i].inputs.logbeta[j]) + fabs (data[i].desired_output.new.logbeta[j])) ||
                  ( fprintf (stderr, "%u %u: %g ?= %g\n", i, j, data[i].inputs.logbeta[j], data[i].desired_output.new.logbeta[j]),
                    0 ));
        }

      for (unsigned int j = 0; j < 5*5; ++j)
        {
          assert (fabs (data[i].inputs.gamma[j] - gammacopy[j]) <=
                  1e-8 * (1 + fabs (data[i].inputs.gamma[j]) + fabs (gammacopy[j])) ||
                  ( fprintf (stderr, "%u %u: %g ?= %g\n", i, j, data[i].inputs.gamma[j], gammacopy[j]),
                    0 ));
        }
    }
}

int 
main (void)
{
  srand48 (69);

  test_to_from_inverse ();
  test_qfunc_f ();
  test_qfunc_df ();
  test_mstep_no_data ();
  test_mstep_symmetric_no_data ();
  test_mstep ();
  test_mstep_no_hyperprior ();

  return 0;
}
